felm.summary.table = function(summaries = list(), 
                              output = "default",
                              fmt = 5,
                              statistic = "std.error",
                              statistic_override = NULL,
                              statistic_vertical = TRUE,
                              conf_level = 0.95,
                              stars = TRUE,
                              coef_map = NULL,
                              coef_omit = NULL,
                              coef_rename = NULL,
                              gof_map = NULL,
                              gof_omit = NULL,
                              add_rows = NULL,
                              title = NULL,
                              notes = NULL,
                              align = NULL,
                              estimate = "estimate",
                              model.names = NULL){
  
  require(lfe)
  require(stringr)
  n.models = length(summaries)
  
  l=lapply(1:n.models, FUN = function(i){
    s = summaries[[i]]
    
    
    mock = data.table(Y = rbinom(s$N,1,.5),
                      T = rbinom(s$N,1,.5))
    
    for(j in 2:nrow(s$coefficients)){
      mock[[paste0('X',j-1)]]=rbinom(s$N,1,.5)
    }
    for(j in 1:length(s$fe)){
      mock[[paste0('FE',j)]]=rbinom(s$N,1,.5)
    }
    
    form = 'Y ~ T'
    if(length(names(mock)[grepl('X',names(mock))])>0){
      covars = paste(c('',names(mock)[grepl('X',names(mock))]), collapse = ' + ')
    }
    if(length(names(mock)[grepl('FE',names(mock))>0])){
      fe = paste(names(mock)[grepl('FE',names(mock))], collapse = ' + ')
      fe = paste('|', fe, '| 0 | 0')
    }
    
    form = paste0(form,covars,fe)
    form=as.formula(form) 
    
    m = felm(form, data = mock)
    
    m$beta[,1]=s$coefficients[,1]
    rownames(m$beta) = rownames(s$coefficients)
    m$se = s$coefficients[,2]
    m$pval = s$coefficients[,4]
    m$N = s$N
    m$residuals[,1] = s$residuals
    m$TSS[1] = sum(s$residuals^2)/(1-s$r.squared)
    m$df = s$df[1]
    m$df.residual = s$df[2]
    m$fe = s$fe
    m$terms = s$terms
    m$hasicpt=s$hasicpt
    
    m$p =  s$p
    
    return(m)
    
    
  })
  
  if(!is.null(model.names)){
    
    names(l) = model.names
  }
  
  out = modelsummary(l,
                     output = output,
                     fmt = fmt,
                     statistic = statistic,
                     statistic_override = statistic_override,
                     statistic_vertical = statistic_vertical,
                     conf_level = conf_level,
                     stars = stars,
                     coef_map = coef_map,
                     coef_omit = coef_omit,
                     coef_rename = coef_rename,
                     gof_map = gof_map,
                     gof_omit = gof_omit,
                     add_rows = add_rows,
                     title = title,
                     notes = notes,
                     align = align,
                     estimate = estimate)
  
  return(out)
  
  
}

#felm.summary.table(summaries = list(ModelDemSpExpDems, ModelDemSpExpReps))
